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Non-Abelian plasma instabilities have been proposed as a possible explanation 
for fast isotropization of the quark-gluon plasma produced in relativistic heavy-ion 
collisions. We study the real-time evolution of these instabilities in non-Abelian 
, plasmas with a momentum-space anisotropy using a hard-loop effective theory that 

I/"") \ is discretized in the velocities of hard particles. We extend our previous results on 

the evolution of the most unstable modes, which are constant in directions transverse 

ff~j ■ to the direction of anisotropy, from gauge group SU(2) to SU(3). We also present 

first full 3+1-dimensional simulation results based on velocity-discretized hard loops. 

\q , In contrast to the effectively 1+1-dimensional transversely constant modes we find 

CN I sub exponential behavior at late times. 

o : 

in ■ 

I. INTRODUCTION 

^ : 

Ph! If a quark-gluon plasma has been produced in the current heavy-ion-collider experiments 

at RHIC, the early onset of hydrodynamic behavior indicates so fast thermalization, or at 
least isotropization, that perturbative scattering processes alone do not seem to be capa- 
ble of explaining it 0, HI IE IH El • A number of nonperturbative scenarios have been put 
. £h ! forward to account for the seemingly inevitable intrinsic strong-coupling nature of the quark- 

gluon plasma 0, 0, H, El • However, a possible nonperturbative effect that does not rely on 
specifically strong-coupling physics, and that has not yet been fully taken into account in 
the available weak-coupling analyses, is the important collective dynamics arising from non- 
Abelian plasma instabilities, which are generalizations of the so-called Weibel or filamentary 
instabilities in ordinary electromagnetic plasmas [10]. Such instabilities are present already 
in collisionless plasmas with any amount of momentum anisotropies [ill . 12l |. and, indeed, 



they have been found to play an important role in the fast isotropization of electromag- 
netic plasmas 0]. Their non-Abelian versions have been pro pos ed to be of relevance for 
the quark-gluon plasma early on by Mrowczyriski and others jlil . EiL Hil 17 , [3, [l^| , and 
specifically as explanation for the recent puzzle of fast apparent thermalization by Arnold, 

Lenaghan, Moore, and Yaffe mini. 



At weak coupling, the asymptotic regime where perturbative QCD becomes justifiable, 
there is a separation of scales between the scale pertaining to hard particles forming the 
nearly collisionless constituents of a plasma and the scale which is down by one power of 
the QCD coupling g, determining for example electrical screening masses. The wave vectors 
and growth rates of Weibel instabilities are of the same order, which is parametrically much 
larger than perturbative collision rates or even the parametrically larger color relaxation 
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rate. The degrees of freedom at the scale gptard are to leading order described by a non- 
local and non-linear effective action called the hard-thermal-loop effective action |22j in the 
case of thermal equilibrium, and this has a simple generalization to plasmas with arbitrary 
momentum-space anisotropy [23l |. 

In the static limit, the anisotropic hard- loop effective action has a potential which is un- 
bounded from below. Using a static approximation and restricting to modes with momenta 
in the direction of anisotropy, upon which the hard-loop effective action becomes Gaussian 
and local, the evolution of plasma instabilities has been studied numerically by Arnold and 
Lenaghan in Ref. 23. Based on these results, these authors have conjectured that non- 



Abelian plasma instabilities Abelianize in the regime of nonperturbatively large amplitudes, 
and thus grow exponentially like their Abelian counterparts as long as the assumptions of the 
hard-loop approximation hold. At parametrically larger gauge-field amplitudes A ~ Phard/g, 
when they begin to have non-perturbatively large effects on the hard particles' trajectories, 
these modes are expected to eventually saturate by isotropizing the hard particle distribu- 
tion. 1 

The full hard-loop effective action, however, is nonlocal and nonlinear. Using an auxiliary- 
field formalism to make it local and discretizing with respect to the velocity space of the 
hard particles, we have carried out a full hard- loop simulation for an SU(2) plasma in 
Ref. |25[ using initial conditions that like those considered in Ref. |24] are constant with 
respect to the spatial directions transverse to the anisotropy axis. These configurations 
evolve 1+1-dimensionally and contain the most unstable modes for a given longitudinal 
momentum. We found that there is subexponential behavior when the instabilities reach 
nonlinear amplitudes, but exponential growth roughly equal to that of the linear regime was 
restored deeper into the nonlinear regime. However, the instabilities of the full hard- loop 
effective theory do not exhibit global Abelianization as in the simplified model of Ref. [24] , 
but only in limited domains with extent comparable to the wavelength of maximal growth. 

In this paper we give a fuller account of these numerical simulations and we present 
their extension to the gauge group SU(3) of QCD. After introducing the local auxiliary- 
field formulation of the hard-loop effective theory in Sect. |HJ we introduce our method 
to discretize hard- loop momenta in Sect. 11111 and compare with continuum results for the 
dispersion laws in Sect. I1VI In Sect. |V1 we study the effectively 1+1-dimensional evolution 
of initial configurations with random fluctuations in the spatial direction of anisotropy, but 
which are constant with respect to transverse directions. We analyse discretization effects 
and the effect of upgrading the gauge group SU(2) to SU(3). In Sect.|Vl]we finally consider 
initial conditions that are non-constant in the transverse direction and thus require full 
3+1-dimensional simulations. While at the beginning of the nonlinear regime we obtain 
results similar to the effectively 1+1-dimensional simulations, we observe subexponential 
behavior deeper into the nonlinear regime, in qualitative agreement with a most recent 3+1- 
dimensional study by Arnold et al. |26( using a different discretization method. Sect. IV 111 
contains our conclusions. 



This important aspect has recently been studied in Ref. |32| f° r non- Abelian plasma instabilities. Our 
objective here is to investigate the hard-loop dynamics preceding (in a weak-coupling regime) the regime 
where backreaction needs to be included. 
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II. HARD-LOOP EFFECTIVE FIELD EQUATIONS FOR ANISOTROPIC 

PLASMAS 

At weak gauge coupling g, there is a separation of scales in hard momenta |p| = p° 
of (ultrarelativistic) plasma constituents, and soft momenta ~ g\p\ pertaining to collective 
dynamics such as Debye screening, finite plasma frequency, and, in anisotropic plasmas, 
magnetic instabilities. The effective field theory for the soft modes that is generated by 
integrating out the hard plasma modes at one-loop order and in the approximation that the 
amplitudes of the soft gauge fields obey -C \p\/g is that of gauge-covariant collisionless 



Boltzmann-Vlasov equations [27[ . In equilibrium, the corresponding (nonlocal) effective ac- 
tion is the so-called hard-thermal-loop effective action E3| which has a simple generalization 
to plasmas with anisotropic momentum distributions [23j. At the expense of introducing a 
continuous set of auxiliary fields, the hard-loop effective field equations can be made local in 
space and time. These field equations then involve an induced current of the form 3 HI 



with a set of auxiliary fields Wp{x\ v) for each spatial unit vector appearing in the velocity 
v ij. = p/y|p| = (l ? v) of a hard (ultrarelativistic) particle with momentum The W's 
satisfy 

[v D{A)]W p {x-v) = F Pl {A)v^ . (2) 

with Dp = 3^ — ig[A^ ■]. 

The equations are closed by the non-Abelian Maxwell equations 

D,(A)F^=f, (3) 

where v = is the Gauss law constraint. 2 

In the following we shall consider the special case of cylindrically symmetric anisotropies 
in momentum space, so that there is only one direction of anisotropy, which in heavy-ion 
collision would be given by the collision axis. We can then parametrize /(p) = /(|p|,p 2 ) 
and write 



df(p) 



df(\p\ 


,p z )p 




d\p 




P 





H 7^ i o. 



dpP <9|p| |p| dp 



zf3 



= -A(IpU z )^-/ 2 (|p|,p 2 )V (4) 

Eq. (j2J) allows us to impose p^Wp = 0, which leads to 

f{x) = \g 2 J -j^-y [~hW\x- v) + f 2 W z {x- v)]. (5) 

Notice that in the isotropic case one has f 2 = 0, and only W° appears, whose equation of 
motion (J2J involves only electric fields. In the anisotropic case, however, W z enters, whose 
equation of motion contains the z component of the Lorentz force. 



Our metric convention is (- 
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The induced current (JIJ is covariantly conserved for arbitrary momentum-&pace distribu- 
tion functions / as can be verified by partial integration with respect to p@ 23] . In the case 
of a parity-invariant function /, /(p) = /(— p), to which we shall restrict ourselves in the 
following, no partial integration is needed and, moreover, the two terms in the current (J3J) 
are covariantly conserved individually: Using (J2J) we have 

D ■ j oc J -^-(f^v, + f 2 F"v y ) = (6) 

where parity invariance of fx makes the v in the first term integrate to zero. In the second 
term f'2 is odd with respect to the ^-direction, so that only the term involving F zz v z could 
contribute, but F zz = 0. 

The dynamical system described by Eqs. (HJ-© has constant total energy of the form 

£ = J d 3 x tr(E 2 + B 2 )\ t + £dt' J d 3 x 2 tr j t , -E t ,. (7) 

The part containing the induced current involves 

4lwy^ w ' {v - D)w - < 8) 

which shows that in the isotropic case (J2 = 0) there is a local, positive definite energy 
contribution from the plasma [29]. However, in the general anisotropic case, positivity is 
lost. This reflects the possibility of plasma instabilities, where energy may be extracted 
from hard particles and deposited into the soft collective fields without bound as long as the 
hard- loop approximation <C \p\/g remains valid. 



III. DISCRETIZED HARD LOOPS 

The structure of Eqs. (HJ-© is such that only W and W z participate nontrivially in the 
dynamical evolution. A closed set of gauge-covariant equations is thus obtained when the 
hard-loop integral in (jSJ is discretized with respect to directions v 

fid = 9 2 J ^^v»[hW»{x) + f 2 W*{x)\ 

^^T,v^[a v W^x) + b v W:(x)], (9) 

where the Af unit vectors v define a partition of the unit sphere in patches of equal area. 
The coefficients a v , b v are constants defined by 

« v = -^§^A(ipUpin (10) 

r°° p 2 dp 

J2^y 



~9 2 I ""£%/ 2 (IpUpI**)- (11) 



5 



Isotropic distribution functions / are characterized by a v 's which are independent of v, and 
vanishing 6 v 's. 

The special choice of Ref. 3] for an anisotropic distribution function, 

f(\p\,p z ) = N(Of iso (p 2 + & 2 z ), (12) 

gives 



(1 + ^ 2)2 , b v = £v z a v , (13) 

where is the Debye mass of the isotropic case £ = 0. If one requires that the number 
density of hard particles is the same for different values £, the normalization factor here is 
= \/l + £• m t ne following we shall often use the abbreviation 

m 2 = iV(0m 2 D , (14) 

leaving open the possibility for different normalizations. 3 

Discretization of the directions v can spoil the reflection properties of f\ t 2 and thus 
automatic covariant conservation of j by virtue of Eq. (JHjl. Covariant current conservation 
is however guaranteed when a v , b v satisfy 

$> v v = 0, £ & v = 0, ]>> v v ± = 0, (15) 

V V V 

where vj_ = v — v z e z , and so we shall restrict ourselves to discretizations of this kind. 

A suitable discretization of the sphere is given by a set of unit vectors v following from 
regular spacing in z and (p according to 

z . = _i + (2z - i)/N z , i = l...N z , i Pj = 27rj/N (fi , j = l...N v . (16) 

The resulting "disco balls" are such that they are covered with M = N z x N v (curved) tiles 
of equal area. For low values of A/", we shall also consider discretized v's pointing to the 
centers of the faces of one of the regular polyhedra with M = 6, 8, 12 or 20. 

Given any set of v, a v , b v , satisfying (|T5|) . one has then to solve (JO]) together with the 
Yang-Mills field equations and 

[v ■ D{A)]WJ{x) = F Pl {A)v\ (3 = 0, 3. (17) 

In temporal axial gauge A = 0, the dynamical variables are Ai, E{ = —Ai, Wq and W£. 
Eq. (|T7j) becomes 

d t WJ{x) = -vT>WJ + F 01 {A)v\ (18) 

The gauge choice A = is in fact merely a convenience for the following numerical simula- 
tions. All equations are gauge covariant, the current is covariantly conserved, and we shall 
also consider exclusively gauge invariant observables in our numerical simulations. 

Through Eq. (JUJ) only a linear combination of W® and W£ participates actively in the 
dynamical evolution. So in addition to the gauge fields we need to consider only the J\f 
auxiliary fields 

W v (x) = a v W°(x) + b v W*(x). (19) 



3 Note that Ref. [ill has AT(£) = 1. 
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The full hard-loop dynamics is then approximated by the following set of matrix-valued 
equations, 



[v ■ D(A)}W V = (a v F 0fl + & V F 2 > M ( 20 ) 
1 

Af 



DM)*"" = f = T?E^Wv, (21) 



which can be systematically improved by increasing Af. 
The Gauss law constraint is explicitly 

Di(A)F i0 = -A(A)i J = j° = ± £ W v . (22) 

In the case of an isotropic plasma, a similar discretization of the hard-loop dynamics which 
uses a discrete set of vectors v has been considered before by Rajantie and Hindmarsh |3(|. A 
different possibility for discretization that has been em ployed p reviously is a decomposition 
of the auxiliary fields W li (x; v) into spherical harmonics [3(J,|31[ and truncating at a maximal 
angular momentum l max . This method has now also been used in the most recent numerical 
study of non-Abelian plasma instabilities by Arnold et al. [26]. Discretizing v is slightly 
simpler, but also more flexible in that it allows one to e.g. improve approximations in highly 
anisotropic cases with cylindrical symmetry by selectively increasing the resolution in the z 
direction more than in the (p direction. 

When the momentum space distribution of hard modes has an oblate form, £ > in 
Eq. (I12D . the momenta of the unstable collective modes form prolate regions with \k z \ > 



|21|. For a given k z , the growth rate is largest for k± = 0, i.e. the constant modes with 
respect to the directions transverse to the direction of momentum space anisotropy. By 
considering initial conditions which are constant in the transverse direction, we can study 
the complete dynamics of these particular unstable modes in a 1+1- dimensional setting 
where all fields are dimensionally reduced, (x) — > (t, z). Then only A z (t, z) plays the role of 
a gauge field, and A X)y (t, z) behave as adjoint matter. 

Using temporal axial gauge, the equations of motion for the dynamical fields can then be 
simplified according to 

d 2 t A x = D 2 z A x -g 2 [A y ,[A y ,A x }}+f 
d 2 t A y = D 2 z A y -g 2 [A x ,[A x ,A y }}+f 
d 2 A z = -ig[A x ,D z A x }-ig[A y ,D z A y }+j z 
(d t + v D)W V = Ov(-v ■ d t A) + K[-d t A z + D z (v x A x + v y A y )], (23) 

where D = V + ig[A,-], V = (0,0,9,), and j = ^vW v . 

In Sect. IVll we shall also consider initial conditions which are not constant with respect 
to the transverse directions. This will then require the full 3+1-dimensional content of 
Eqs. (120]) and flUJ). 



IV. COMPARISON OF DISPERSION LAWS 

Before turning to the numerical evaluation of the above equations by also discretizing 
space and time, we shall study the dispersion laws of linearized modes, and the consequences 
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of replacing the continuous set of fields W^(x, v) by a finite number W@(x). Because we shall 
be interested mainly in the unstable modes, we shall concentrate on wave-vectors parallel 
to the anisotropy axis, which includes the most unstable modes. With k % = k8 l t here are 
then only two possible modes for gauge fields, transverse and longitudinal ones 11. 12l . 



A. Continuum results 



For anisotropic distribution functions ()12|) . analytical results for the gluonic self-energies 
have been obtained in Ref. 11, which we recapitulate for the special case of longitudinal 



wave vectors. 



1. Transverse polarizations 



With rj = u/k, the transverse gluon self-energy reads [1 



n, 



2 

m 



1 + V 2 + £(-1 + (6 + OV 2 - (1 - OV 4 )) arctan J£ 



rj — 1 + ie / 

where we recall that m 2 = N(^)m 2 D (C, =0). 



(24) 
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The limit k — > or 11 — > 00 determines the transverse plasma frequency as 

w it = ^(y 1 + arctan V^J m2 - ( 25 ) 

For large momenta, the transverse excitations tend to the asymptotic mass given by 

mle/m 2 = — -= arctan J~£. (26) 
2v4 

The static limit on the other hand gives a mass squared which is determined by the 
relation 

< = -V 2 = ~Hi,t ■ (27) 
For negative £, this indicates screening of magnetostatic fields, while for positive £ there is 
an instability for all momenta k < ix. 

2. Longitudinal polarizations 
In the longitudinal sector we have [l2T | 

1 +0C 1 - £r/ 2 ) arctan ^ 



2 2 
77 m 

il = 



+V / ?((l + ^ 2 )-(l + 0^n r? + l^ e 



(28) 



4 



For — 1 < £ < one has to replace arctan v 7 ^ by — ^=artanbV~? 
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The plasma frequency in the longitudinal sector is generally different from the one in the 
transverse sector when £ 7^ and reads 



^lu = 7^ {—Jf arctan \/£ - 1 ) m2 • ( 29 ) 



The static limit of He/t] gives the electric screening mass squared 
2 1 1 1 s „ /7 1 1 I _2 _2 1 ,2 



m el = - ^— ^arctan ^ + ll m = m + f . (30) 

This parameter does not change sign for all —1 < £ < 00, so in the present case of wave 
vectors parallel to the anisotropy direction all instabilities are purely transverse. 



B. v-discretized hard loops 

After discretizing the space of velocities and approximating the fields W@(x, v) by a finite 
set of field W£(x), we search for solutions of Eqs. (|23|) where all fields point in a fixed color 
direction and which correspond to a single Fourier mode in z-direction: 

A% z) = e 7 ¥ fc V, W°' z (t, z) = e^e ikz w° v z . (31) 

Here k will usually be a real number, while 7 can be either real or imaginary, corresponding 
to exponential or oscillatory behavior, respectively. The dynamics is then effectively Abelian 
and linearized, and Eqs. (j2HJ) reduce to 

7 y + k\\ = ^7 E + 6v<) (32) 

(y + ikv g )wl = -7«V (33) 
+ = - 7 e z + ifo/ei (34) 

with = (e x , e y , 0). The Gauss law constraint reads 

- zA^e* = E( a v^v + M£)- (35) 

Next we specialize to the anisotropic distribution function (|12j). for which the coefficients 
a v , b v are given by (fT3j) . and compare with the above continuum results. With discretized 
v's we obtain algebraic equations that can be solved in closed form for the regular polyhedra 
with M = 6, 8, 12, 20. Taking the v's to point to the centers of the faces of these polyhedra, 
we consider two different orientations for each polyhedron: one where the discrete rotational 
symmetry around the z-axis is maximal, and one where this symmetry is only Z3. The two 
choices sample in particular different sets of ^-components of hard momenta (in Table |H] in 
Appendix IA 21 the number of distinct ^-components is indicated by N z ). We shall also work 
out the dispersion laws in a few cases of larger Af using Eqs. (|lb|). 
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1. Transverse polarizations 

Equations give algebraic relations between 7 and k, while the Gauss law constraint 
is automatically fulfilled for reflection-symmetric distribution functions, which lead to the 
restrictions (fT5J) . 

For the cube (Af = 6) and Z4 symmetry about the z-axis, the solution is simply 

7 2 = -u 2 = -(k 2 + m 2 /3), (36) 

corresponding to a propagating mode with momentum-independent mass equal to the 
isotropic plasma frequency In this case, the ^-dependence drops out completely (except 
for a possible normalization of m according to ([14))) . The result (|3*Sjl is in fact identical to 
the one obtained by discretization through spherical harmonics and truncating at angular 
momentum l max = 1 |3~D ]. 

For Af = 6 and Z 3 symmetry about the z-axis, the solution is more complicated, and 
happens to coincide with the Af = 8, Z\ case. These two cases give the biquadratic equation 

2 2 _ m 2 - 7 2 + £fc 2 /3 

t + k ~ (TTW 3 7 2 + ■ {67) 



This equation obviously has solutions for real 7 (i.e. instabilities) only when £ > 0, as is the 
case for the continuum solutions. There are two branches of solutions given by 



2 2 1 

7 =-u = -- 

o 



a + Ak 2 ± J a 2 + (8 + 4f)afc 2 + 4fc 4 



(38) 



The upper sign corresponds to propagating transverse plasmons. As can be seen from Table 
ITTl in Appendix IA 21 there is considerable improvement compared to the simple dispersion 
law ()36|) of the M — 6, Z± case — the asymptotic mass and its small-£ correction are even 
exact. The lower sign in (|38j) contains the instability for £ > 0, and oscillatory behavior for 
k > n. Fig. shows a plot of the solutions of (|3*5|l and compares it with the continuum 
result following from ([24)1 . 

The M = 8, Z 3 case and the N = 12, Z 5 case also lead to biquadratic equations, each 
with somewhat different coefficients, and the small-£ behavior can be read from Table ITTl 

The Af = 12, Z 3 case and the icosahedron, Af = 20, on the other hand, lead to bicubic 
equations, with Af = 12, Z 3 and Af = 20, Z 5 sharing the same equation 5 , namely 

2 2 2 (45 7 4 ~ A: 4 Q(45 + 42£ + 13£ 2 ) + 3 7 2 fc 2 (315 + 255£ + 93£ 2 - 7£ 3 ) 

7 (45 + 30£ + £ 2 ) 2 (457 4 + 307 2 P + fc 4 ) ' 1 j 

The case Af = 20, Z 3 gives a similar equation but with different coefficients. The solutions 
to these equations can be given in closed form, but are too unwieldy to be reproduced here. 
There are now three distinct solutions one of which corresponds to propagating transverse 
plasmons. 



5 As is well known, the regular polyhedra Af = 6, 8 and Af — 12, 20 form dual pairs, however this does not 
seem to explain the coincidence of the dispersion laws of Af = 6, Z3 with Af = 8, Z4 and of Af = 12, Z3 
with Af — 20, Z5, since duality would map orientations with rotational symmetry Z3 to such with Z3. 
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(a) N = 8 (b) N = 20 

FIG. 1: The transverse dispersion laws u(k) for anisotropy parameter £ = 1 in the continuum 
(full lines) compared with those of v-discretized hard loops for the octahedron (M = 8, Z4) and 
the icosahedron (N = 20, Z5) (dashed lines). The unstable modes appear at small space-like 
momenta where the plot is in terms of j(k) instead of iv(k). The additional space-like modes are 
approximations of the Landau cut by simple poles whose number increases with N . 



Except for the cube with Z 4 symmetry, Eq. (JHrjj) . there is always exactly one solution 
which contains an instability for £ > and k < fi. A significant difference to the continuum 
result is however the absence of Landau damping. Instead, there are additional modes for 
spacelike momenta, and their number increases with A/", approximating the logarithmic cut 
by more and more poles at spacelike momenta. This is analogous to what has been observed 
in Ref. [3]J for discretizations using spherical harmonics. However, in the anisotropic case, 
one of these spacelike modes becomes unstable for £ > 0. In the continuum limit, the 
unstable part (k < fi) remains on the physical sheet, whereas its stable, oscillatory part 
moves to the unphysical sheet, cf. Ref. 



In Fig. we show the transverse dispersion laws for £ = 1 following from (JH9"j) . which 
has two spacelike branches. The comparison with the continuum case shows that both the 
timelike propagating modes and the instabilities are quite well reproduced quantitatively. 

Our main interest will clearly be the plasma instabilities. Since all the actual and addi- 
tional stable modes are just oscillatory, we do not expect them to play an important role in 
the evolution of the instabilities. For us the most important aspect of the dispersion laws 
will be the magnitude of 7(fc) for the unstable modes k < fi. The continuum result for 
'y(k) has the property that it vanishes both at k = and at k = //, in contrast with the 
toy model studied in Ref. [2^. This feature is correctly reproduced in our discretization of 
the hard-loop dynamics of anisotropic plasmas. In Figs. |2]we compare ^y(k) of the various 
polyhedra (with maximal Z v ) with the exact continuum result for some values of £. From 
this we expect that all the polyhedral approximations should be qualitatively correct, with 
Af = 20 being a quantitatively good approximation for £ < 10. 

For higher Af we use the "disco-ball" discretization . Fig. |3] compares the growth 
rates j(k) for full and discretized hard loops with Af = Af z x Af v = 20 x 5 = 100 and 
£ = 10, which shows that the latter give very accurate approximations for Af ^ 100. Also 
the dispersion laws of propagating stable modes are accurately reproduced. For example, 
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l(k)lm i=\ i=\q ^=100 

o. I . , , , , , — . . , , , , , . 




0.1 0.2 0.3 0.4 0.5 0.2 0.4 0.6 0.8 1 0.25 0.5 0.75 1 1.25 1.5 1.75 



klm klm klm 

FIG. 2: The growth rate j(k) of unstable modes for the polyhedral approximations N = 8, 12, 20 
and maximal £L, (dashed lines with smaller and smaller dashes) in comparison with the continuum 
result (black lines) for £ = 1, 10, and 100. 

0.12 

o.io 

-S 0.08 
^ 0.06 

0.04 

0.02 

0.2 0.4 , , 0.6 0.8 1 

klm 

FIG. 3: The growth rate j(k) of unstable modes for the discretized hard loops with M = 100 
(dashed line) in comparison with the continuum result (solid line) for £ = 10. The precise values 
of wave number k* and growth rate 7* corresponding to the mode of maximal growth can be read 
from Table [H 

with £ = 10, jV = 100, the continuum result for the asymptotic mass (|A3|) of the transverse 
modes with |k| ^> /x is reproduced with an error of only 0.017% 

In the following numerical simulations, we shall employ the disco-ball discretization with 
M > 100, choosing N z 3> N v to take into account that the distribution function varies more 
rapidly in the z than in the 92-direction. When considering smaller Af, we shall instead 
take regular polyhedra, which as discussed further in Appendix IA 2|. seem to be superior 
approximations for M < 20. Before turning to these numerical simulations, we complete 
this section by briefly considering also the (stable) longitudinal modes. 




2. Longitudinal polarizations 

For longitudinal polarizations (e 2 7^ 0) the Gauss law constraint is nontrivial, but it 
leads to the same equation as the linearized field equation (J32~j) . For the regular polyhedra, 
the resulting algebraic relations between 7 and k can be solved in closed form, and one finds 
that as expected there are only stable solutions — 7 2 = uj 2 > 0. 
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(a) N = 8 



klm 

(b) M = 20 



FIG. 4: The longitudinal dispersion laws uo{k) for £ = 1 in the continuum (full lines) and for the 
octahedron {N = 8, Z4) and the icosahedron (AT = 20, Z5) (dashed lines). 

For M = 6 and Z4 symmetry about the z-axis, the solution is 



u 2 = 



m 



3(1 + 



+ k 2 , 



(40) 



corresponding to a propagating mode with momentum-independent mass which at £ = 
equals the correct plasma frequency. In the static limit, the Debye mass is off by a factor 
1/V3. 

The cases M = 6 with Z3 symmetry about the z-axis and M = 8 with Z4 symmetry have 
a similarly simple solution, which however reproduces both the plasma frequency and the 
Debye mass in the isotropic limit through the dispersion law 



2 _ m 2 (l + Q k 2 
U 3(l + £/3) 2+ 3- 



(41) 



At larger momenta the solution becomes spacelike, whereas the continuum dispersion law 
would approach the light-cone exponentially in the timelike domain. 

With increasing J\f, the continuum result for the longitudinal dispersion relations is ap- 
proximated from below (i.e. spacelike momenta) by the mode with largest frequency. The 
Landau cut at spacelike momenta is approximated by an increasing number of poles, one of 
which is connected to the time-like plasmon mode. 6 For M = 12 and 20, there are again two 
spacelike branches, now given by biquadratic equations. The cases M = 20, Z 5 and M = 12, 
Z 3 again share the same equation, which reads 



45cu 4 - 30u 2 k 2 + k = 



135m 2 (l + 
45 + 30£ + £ 2 



[5co 2 (45 + 6£ + e) - k 2 (15 + lOf + 3£ 2 )]. (42) 



6 In the continuum, the longitudinal plasmons have effectively a finite range of possible momenta because 
the residue vanishes exponentially as the light-cone is approached. 
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V. 1+1-DIMENSIONAL LATTICE SIMULATIONS 

In the previous section we have seen that for oblate momentum distributions of plasma 
constituents, £ > 0, there are unstable modes for \k z \ < /i and sufficiently small k±. In 
the linear regime, where field amplitudes A are much smaller than k/g, time evolution is 
determined by the above dispersion laws and all modes evolve independently. Unstable 
modes grow exponentially with growth rate 7(k), which for a given value of k z is largest 
when k± = 0. We denote the scale of maximal growth by k* and the corresponding maximal 
growth rate by 7*. This scale is of the same order of magnitude as the Debye screening mass 
or the thermal mass of propagating plasmons, namely ~ gphard, which at weak coupling is 
much larger than the color relaxation scale ~ g 2 phard or the scale of large-angle scattering 
rates ~ g 4 Phard- This is what makes plasma instabilities a prime candidate for the mechanism 
of fast isotropization in a weakly coupled plasma. 

In the following we shall study the evolution of collective soft fields starting from tiny 
random fluctuations. If the initial field amplitudes are sufficiently small to give a large 
number of e-folds in the linear regime, this will lead to field configurations that are dominated 
by the modes of largest growth. For £ > this will favor modes with \k z \ ~ fc* and k± « 0. 
The special initial conditions with strictly k± = 0, i.e. constant modes with respect to the 
spatial directions transverse to the anisotropy direction should therefore be an idealization 
of particular interest. All fields are then functions of only one spatial coordinate, and the 
equations of motion are simplified according to Eqs. (j^j) . 

In Ref. |24|, non-Abelian dynamics with such initial conditions was studied numerically 
with a drastically simplified model for the induced current (jl]). namely 

= n 2 A aj a = x,y. (43) 

As was shown in Ref. jH], this correctly reflects the static limit of the anisotropic hard- 
loop effective action for fields that vary only in the anisotropy direction, but it neglects its 
general frequency dependence and dynamical nonlinearity. Already at the linearized level, 
the former complication means that modes with vanishing wave vector, k = 0, are stable, 
see Figs. El and El whereas the toy model (}4*3*j) implies a growth rate 7 = vV 2 — k 2 , which 
is maximal at k = 0, where it equals /z. As Figs. E] and HU show, the anisotropic distribution 
function (fl"2"|) with £ > leads to 7* < k* < /i instead. 

In simulations of the toy model (|4*H|) . the authors of Ref. [24j have found that unstable non- 
Abelian modes, which are constant modes with respect to transverse space, might behave 
very similarly to Abelian Weibel instabilities also in the nonlinear regime. In the latter, they 
observed rapid Abelianization, both locally and globally. 

Using the v-discretized equations of motion, (|2Uj) and (|2ip. we have extended the lattice 



simulations of Ref. 24] to the full hard-loop effective theory, with the main findings presented 
in Ref. [2^. In the 1+1- dimensional situation we found some modification of the evolution 
of the instabilities when they first reach nonlinear size. Subsequently, growth with rates 
close to the Abelian growth rate was restored, however. 

In the most recent paper [2^ a similar analysis was performed, using a different dis- 
cretization method for the hard-loop effective theory. While confirming our 1+1-dimensional 
findings, they also considered random initial fluctuations in all spatial directions and the 
resulting 3+1-dimensional evolution with the conclusion that the late-time behavior, deeper 
in the nonlinear regime, is in fact modified importantly. In this section we shall concen- 
trate on the 1+1-dimensional situation, giving a detailed account of the results presented in 
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TABLE I: Debye mass m e i, magnetic instability scale [a, wave vector fc* and growth rate 7* for 
the modes of maximal growth, transverse and longitudinal plasma frequency, and the asymptotic 
mass of transverse gluons for the isotropic case and the anisotropic case assumed in the numerical 
simulations. 





m e i/m 


fj,/m 




7*/m 


u p i,t/m 


Upi,e/m 


rrioo/m 





1 








0.577350 


0.577350 


0.707107 


10 


1.64296 


1.07225 


0.398025 


0.120289 


0.339075 


0.412228 


0.447144 



Ref. [25]. After extending this also to gauge group SU(3), we shall give our first results for 
a 3+1-dimensional simulation of the v-discretized hard- loop effective equations in Sect. IV11 

Like Ref. |2J] we work in temporal axial gauge, A = 0, and take initial conditions cor- 
responding to small random chromoelectric fields E = —d t A with polarization transverse 
to the z-axis, and all other fields vanishing, which in our case includes the auxiliary fields 
W v - This initial condition satisfies the Gauss law constraint, D ■ E = A/" -1 J2 V W v , whose 
continued fulfilment is monitored in the simulation, but not enforced. As a further non- 
trivial check of our numerics, we tracked conservation of the total energy (|7J), stopping the 
simulation when this signalled loss of accuracy. The lattice versions of the 1+1 dimensional 
equations of motion that we have employed are given in detail in Appendix [BJ 

In the 1+1-dimensional simulations below we used anisotropy parameter £ = 10. We shall 
express all dimensionful quantities by the asymptotic mass of transverse excitations, whose 
relation to the mass scale m appearing in the coefficients (|13|) is given in Table |1] While the 
coupling constant g could be absorbed by redefinitions of our dynamical quantities, we shall 
make their appearance explicit in our final results. 

The lattice spacing used in the simulations presented in this section was chosen such that 
a = 0.0707771^. We have taken a one-dimensional spatial lattice with periodic boundary 
conditions and 5,000 sites, so that the physical size is L m 350m {X) 1 . We shall show results 
for the icosahedron Af = 20, Z$, and disco-ball discretizations M = 100, 800, and 2000, with 
N z x = 20 x 5, 50 x 16, and 100 x 20, respectively. In the leapfrog algorithm specified in 
Appendix IB1 we used time steps e/a from 1/100 to 1/250. The initial conditions are given in 
Eq. (jB10|) . and the choice of the parameter a therein corresponds to starting with a random 
initial seed chromoelectric field of root-mean-square amplitude 0.012 m^/g. 

In the numerical simulations, the Gauss law constraint, which is satisfied by the initial 
conditions, remained satisfied within machine accuracy throughout. As a further check, we 
monitored the total energy ((7|) and stopped the simulation when this had 0(1) violations. 
For all but the largest times, this quantity was conserved within less than 1%; only at the 
very end when the energy in the various field components had grown by factors larger than 
10 6 was there energy violation at the percentage level, which quickly thereafter exploded, 
signalling uncontrollable discretization errors. 
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A. Currents and energies 

Fig. El shows the evolution of the average root-mean-square transverse and longitudinal 
currents and charge density defined by 

1/2 

(44) 

comparing the cases M = 100, 800, and 2000. The transverse current grows exponentially 
with a growth rate that is most of the time only slightly below 7*, except for a transitory 
reduction at the beginning of the nonlinear regime, when j rms ~ m^/g. Our initial condition 
(only nonvanishing transverse electric fields at t — 0) give rise to longitudinal currents and 
charge densities only through non-Abelian interactions (in the Abelian regime all modes 
would decouple and evolve independently according to the linearized dispersion laws). The 
non-Abelian interactions naturally lead to growth rates for j z and j° which are double that 
of J -1 , and such a behavior is seen in Fig. El after some initial delay caused by the fact that the 
unstable modes first have to grow out of the initial random mixture of stable and unstable 
ones. The different cases jV = 100, 800, and 2000 lead to only slightly different behavior in 
the nonlinear regime, where higher M seems to be required to capture the precise variations 
of the subdominant components j z and j°. 



Jrms 



£ ^2tr((j^' ) 2 ) 




FIG. 5: The average current j^ ms (black), j* ms (blue) and rms charge density j® ms (red) for 
N\\r = M = 100, 800, 2000 using a collection of 4 runs each with different random initial conditions. 



Closely related to this, Fig. |U] shows how the exponentially growing energy that is trans- 
ferred from hard to soft scales gets distributed among chromomagnetic and chromoelectric 
collective fields. The dominant contribution is in transverse magnetic fields, and it grows 
roughly with the maximum rate 7* both in the linear as well as in the highly nonlinear 
regime, with a transitory slowdown in between. Transverse electric fields behave similarly, 
and are suppressed by a factor of the order of (7*/fc*) 2 . Notice that in the model of Ref. 24] 
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the situation would be different: Because k* is zero there, the dominant energy component is 
then from transverse electric fields, whereas the relative importance of magnetic fields drops 
with time. In the 1+1-dimensional evolution, the transverse magnetic field component is 
the dominant one, both in the linear and in the nonlinear regime. 

The appearance of longitudinal contributions, which are absent in the initial conditions 
we have chosen, is again a purely non-Abelian effect. While completely negligible at first, 
they have a growth rate which is double the one in the transverse sector, and they begin 
to catch up with the latter when the nonlinear regime is reached first. At this point, the 
general growth stalls for a time of order 7" 1 . When the transverse magnetic field resumes 
its growth (with the transverse electric field following with some delay), the respective 
longitudinal components drop for some time before also starting to grow again. This dip in 
the subdominant longitudinal components turns out to depend on A/", whereas the overall 
behavior is already reproduced well by the rather small value of M = 20 (the icosahedron) . 

Deeper into the nonlinear regime, the (transversely constant) unstable modes grow at 
roughly equal rates. There the results are only weakly dependent on the degree of hard-loop 
discretization, Af. 



B. Color correlation 

Had we assumed initial conditions where all fields would point in the same color direction 
all over the spatial lattice, we would have found only strictly exponential growth once the 
stable modes have become of negligible importance. The behavior would then be exactly 
the same as with Abelian Weibel instabilities, which grow exponentially until they come 
into conflict with the assumptions of the hard-loop approximation (namely that the fields 
have only small effects on the trajectories of the hard particles), whereupon isotropization 
is supposed to set in. 

Apart from a brief transitory slowdown, the 1+1-dimensional simulations thus evolve 
similarly to Abelian instabilities both in the linear and the strongly nonlinear domain (at 
least as far as the dominant transverse components are concerned). A measure of local 
"non-Abelianness" can be defined by 7 

r<U\ - [ L ^ {tr ((i\j x ,jy]) 2 )} 1/2 a ,s 

C[J] ~ Jo 1 MlTif) ' (45) 

Indeed, when the nonlinear regime is entered, this quantity is found to decay roughly expo- 
nentially (with some hiccup brought out only with sufficiently large A/"), but not as fast as 
it was found in the toy model of Ref. 24] . 



Ref. 24] also observed a concurrent global Abelianization, by comparing the correlation 
among parallel transported spatially separated commutators to the correlation of parallel 
transported fields. Generalizing the quantity used in Ref. [13|, 8 we define 

v m = N "~ l [ L dz trmUz + OMiz + ^z^jz)]) 2 } 

XAU 2N C Jo L tr{jf(z + 0}tr{j?(z)} 1 ° J 



7 This definition coincides with the definition in Ref. 24] in the simplifying case (|43(l , but is gauge invariant 
in the full theory also when the restriction to 1+1 dimensional configurations is removed. 

8 We use currents instead of transverse gauge field components, because they transform like adjoint matter 
when the restriction to 1+1-dimensional situations is removed. 
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FIG. 6: Average energy densities £ in transverse/longitudinal chromomagnetic/electric fields and 
the total energy density contributed by hard particles, £ (HL), for various values of M. In the 
case J\f = 20, where there is a noticeable difference between 7* at the discretized level from 
the continuum value (see middle plot in Fig. [2J), we have rescaled the factor in front of t 
by / = 7*(AA = 20)/7* = 1.29336. The horizontal line terminating in fluctuations is the total 
conserved energy (J7|). 



where U(z', z) is the adjoint-representation parallel transport from z to z' . When colors are 
completely uncorrected over a distance £, this quantity equals unity; if they point in the 



same direction, this quantity vanishes. Following Ref. 24] we define the "Abelianization 



correlation length" £4 as the smallest distance where xa is larger than 1/2, 

^ ] = ^,P- (47) 

This we compare with a general correlation length, not focussing on color, defined through 
the gauge invariant function 



J L f t i{ 3l {z + j)U{z + ^z ) 3l {z)} 
/ L f tr {;,(*) 
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This function now vanishes when fields are uncorrected over a distance £, and it is normal- 
ized such that x(0) = 1. We thus define the general correlation length through 

= min (£). (49) 

The evolution of the general and the Abelianization correlation length is shown in Fig. |HJ 
In contrast to the toy model of Ref. j^, where £a was found to rapidly grow to full lattice 
size when C begins its decay (while no such behavior was seen in general correlation lengths), 
we find that Abelianization takes place only over domains of limited sizes, see Fig. |H| Deeper 
in the nonlinear regime, £4 is seen to have magnitudes comparable to the wavelength of the 
modes of maximal growth, 9 A* = 27r/A;* ~ 7.06 m" 1 . 

In this regime, the different simulations with different M agree fairly well, only in the 
transition region between linear and nonlinear regime it is essential to have higher values of 
M to bring out an initial peak in the evolution of £a- 



C. Extension to gauge group SU(3) 

The above numerical simulations have been performed with gauge group SU(2), which 
is the most economical in computer time. There is no particular reason why non-Abelian 
plasma instabilities should strongly depend on the gauge group (in particular when local 
Abelianization occurs). Nevertheless, this needs to be checked. In this section we discuss 
the changes which occur when the gauge group is taken to be SU(3) instead of SU(2). 
Functionally, the algorithm is exactly the same with the only difference being the underlying 



9 Note again that in the model of Ref. |24| fc* vanishes, which is presumably responsible for the different 
global behavior. 




FIG. 8: The Abelianization correlation length £a[j] compared to the ordinary correlation length 
£[?'] for N = 100,800,2000. 




FIG. 9: The average current j^ s , j^ms an d rms charge density j^ ms with gauge group SU(3) and 
SU(2), for one M = 100 run each. 



matrix representation for the fields. Figures El - ^] show the effect of upgrading the gauge 
group SU(2) to SU(3) for M = 100. 

In Figs. Eland^l we compare the results for the various current and energy components 
in the case of gauge group SU(3) with those for SU(2), both with M = 100. Only minor 
differences are seen, with a little higher energies being reached in the nonlinear regime in 
SU(3), but that would need more extensive tests to be taken as a conclusion. 

The amount of local and global Abelianization is shown in Figs. ^2 and ^1 respectively. 
Again, only small differences can be observed. 
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SU(2) SU(3) 




FIG. 10: SU(3) vs. SU(2) comparison of the average energy densities £ in transverse/longitudinal 
chromomagnetic/electric fields and the total energy density contributed by hard particles, £(HL). 




FIG. 11: SU(3) vs. SU(2) comparison of the local Abelianization measure C of commutators [j x , j y ]. 



VI. 3+1-DIMENSIONAL LATTICE SIMULATIONS 

We finally turn to the question to what extent relaxing strict translation invariance in 
transverse directions, which led to the great simplification of effectively 1+1-dimensional 
dynamics, does modify the evolution of plasma instabilities. In the linear regime, the time 
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FIG. 12: SU(3) vs. SU(2) comparison of the Abelianization correlation length £a[j] together with 
the ordinary correlation length 

evolution favors the modes which are constant with respect to transverse directions, because 
these are the ones of largest growth. Thus, if one starts out with arbitrarily small seed 
fields, transverse translational symmetry will be arbitrarily well prepared for the subsequent 
nonlinear evolution (after correspondingly arbitrarily long Abelian evolution of course). This 
is the case that was covered by the previous section. However, generic initial conditions that 
take a finite time to evolve into the nonlinear regime (or start out there) will always have 
some breaking of transverse translational symmetry, and then one has to simulate full 3+1 
dimensional dynamics to take this into account. 

Figures ITBTfT^l show our first results of such 3+1-dimensional lattice simulations. In these 
simulations we have no longer set up initial conditions corresponding only to transverse 
electric fields, but formulated these initial conditions in terms of the W fields, see Eq. flB20|) . 

We have investigated both asymmetrical and symmetrical lattices, but here we only 
present results from symmetric ones. In Fig. EH we plot the resulting energy densities for 
a 96 3 lattice and the icosahedral approximation, M = 20, compared to a one-dimensional 
simulation using the same parameters. As can be seen from this figure, differences from the 
one-dimensional simulations begin to appear at m^t ~ 35 which is when the field amplitudes 
have reached order m^/g. The chosen lattice spacing corresponds to a = 0.28m^ so that 
the physical volume of the lattice is roughly (28/moo) 3 . In Fig. El we compare the energy 
densities obtained using Af = 20, 100, 200, namely the icosahedron and disco-balls with 
■A4 x J\f<p = 20 x 5 and 25 x 8, with lattice sizes 96 3 , 88 3 , and 69 3 , respectively. The thick 
lines are the energy extracted from the hard particles (HL energy) and the other various 
components of the total energy are shown for comparison but not labelled explicitly. As 
can be seen from this figure there is general agreement between all three discretizations, the 
chief difference being that for late times the then approximately linear growth rate seems 
to decrease with increasing Af. The fluctuations in the overall amplitude are due to the 
different random initial conditions used for each run. 

In the linear regime, i.e. field amplitudes <C Tn^/g, the differences are only due to the 
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FIG. 13: Average energy densities £ in transverse/longitudinal chromomagnetic/electric fields and 
the total energy density contributed by hard particles, £(HL) for 3+1-dimensional J\f = 20. 




100 



m«,t 



FIG. 14: Comparison of average energy densities £ for 3+1-dimensional simulations with J\f = 
20, 100,200 on 96 3 ,88 3 ,69 3 lattices. Thick lines correspond to £(HL) and thin lines are the com- 
ponents which can be inferred from Fig. 1151 In the case N = 20 we have rescaled — ► f m oo 
with / = 7*(AA = 20)/7* = 1.29336 in order to properly normalize the results for comparison. In 
addition, we shifted the J\f = 20 data 0.9 units in m^t to the right and the N = 100 by 0.1 units 
in moot to the left in order to roughly align the exponential-growth phase in the Abelian regime. 
Inset shows late-time behavior of the hard-loop energy density with a linear scale. 



different amount of initial energy densities in longitudinal and transverse components. In 
both the 3+1-dimensional and the 1+1-dimensional case, the dominant component is from 
transverse magnetic fields, followed by transverse electric fields (suppressed by a factor of 
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FIG. 15: Average energy densities £ in transverse/longitudinal chromomagnetic/electric fields and 
the total energy density contributed by hard particles, £(HL) for 3+1-dimensional M = 200, in 
linear scale. 




FIG. 16: The relative rms size of commutators Cfj = tr([jj, jj]) 2 /{trjf trj'J) as a function of time 
compared for 3+1-dimensional versus 1+1-dimensional evolution with M = 200. 

(7*/A;*) 2 < 1. Longitudinal components are initially absent in the 1+1-dimensional prepa- 
ration of the system, but not in the full 3+1-dimensional case because of the different initial 
conditions (|B20|) . Nevertheless, at the later stages of the linear evolution, they behave 
rather similarly. The regime where specifically non-Abelian effects come into the play first 
is also very similar in the two cases, and remain so for a while. However, whereas the 1+1- 
dimensional modes eventually recover their initial growth rates of the linear regime, the 
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transversely nonconstant modes seem to be the determining factor deeper into the nonlinear 
regime, and these behave rather differently from the former. In Figs. EH and one sees 
subexponential behavior at late times, which in fact seems to tend to a linear growth at 
large times. This can be seen most clearly in Fig. where we plot the components of the 
energy density at late times on a linear scale for a 69 3 lattice and M = 200. The chosen 
lattice spacing corresponds to a = 0.28m^ so that the physical volume of the lattice is 
roughly (19/moo) 3 . 

This finding is in agreement with the most recent results by Arnold, Moore, and Yaffe (26^ . 
using a different discretization method for the hard-loop effective theory (namely truncation 
of expansions in spherical harmonics) . While Ref . [26| presented a very extensive analysis of 
possible systematic uncertainties, our results, obtained with cruder discretizations, are still 
somewhat preliminary and require further systematic studies. Nevertheless, we already find 
agreement also in some of the finer details, such as the relative importance of longitudinal 
field components in the late time evolution. 

In Fig. El we moreover display results for the normalized size of the commutators Cfj = 
tr([jj, jj}) 2 / {trjf tijj) which give a measure of local Abelianization with respect to the three 
possible spatial directions. As can be seen from this figure there is initially Abelianization 
in the xy commutator starting at m^t ~ 40, while the other two commutators show no 
evidence for Abelianization. At later times {m^t >~ 60), the xy commutator becomes less 
and less abelian and at late times, in the linear-growth phase, the soft-current commutators 
are nearly isotropic. 

VII. CONCLUSION AND OUTLOOK 

In this paper we have presented detailed results from our numerical studies of discretized 
hard loops applied to the time-evolution of quark-gluon plasma instabilities due to mo- 
mentum space anisotropy. We extended our previous effective 1+1 dimensional analysis to 
include the SU(3) gauge group finding that the results obtained for SU(3) were qualitatively 
similar to those of our SU(2) simulations. In addition, we have presented first results for the 
full 3+1 dimensional evolution of the system using our method to discretize the hard-loop 
effective theory. We tentatively confirm the results of Arnold, Moore, and Yaffe [26| who 
find linear instead of exponential growth of the hard-loop energy at late times. 

The fact that 3+1 dimensional nonabelian plasma instabilities seem to nearly saturate 
in the nonlinear regime where field amplitudes are of order m^j g makes it even more 
interesting and also self-consistent to study their dynamics and their role in isotropization 
and thermalization of quark-gluon plasma within discretized hard-loop approximations. A 
more extensive and thorough study of 3+1 dimensional quark-gluon-plasma instabilities will 
be the subject of a forthcoming publication. 
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APPENDIX A: SMALL £ EXPANSION OF DISPERSION LAWS 



1. Continuum results 



The small-f expansion of the continuum transverse plasma frequency (|25|) reads 

1 2 3 



2 / 2 



— — £ + — £ 2 + Off 3 ) 

3 15 ^ 35 ^ VW 



(Al) 



PL 



l>2 



P.i 



(The coefficients pk are compared with the v-discretized results in Table ITT]) . 
Propagating modes with small k are characterized by 

^(k) = ^ + 0(e)] k 2 + 0(k*). 



6 175 

91 <32 

The small-f expansion of the asymptotic mass (j2*o]) of the transverse excitations reads 

1 1 1 



(A2) 



2/2 



2 6 ^ 



''2 



10 



r-.i 



(A3) 



The plasma frequency in the longitudinal sector, Eq 

uli t e/m 2 



, has the expansion 



3 15 ^ 35 S ^ ' 



(A4) 



*1 «2 



2. Discretized hard loops 



Table HT1 shows how well the coefficients PkiQki^k of the small-f expansion of the above 
parameters in the dispersion laws are reproduced with v-discretized hard loops. In the case 
of regular polyhedra, more and more coefficients Pk,Qk, ^k become exact or come closer to the 
continuum result as M increases. For the two pairs of polyhedra where the dispersion laws 
coincide, there is particular good agreement with the continuum results for the asymptotic 
transverse masses, however at the expense of having one coefficient with a wrong sign in the 
longitudinal sector. 

The lower part of Table ITD performs the same 10 comparison for a few cases of the "disco- 
ball" discretization (|16|) with regular spacing in z and (p. For M < 20, the regular polyhedra 
turn out to superior in the small-f expansion. However, the disco-ball discretization allows 
us to increase M without limit, and while (with the exception of r\) no coefficient becomes 
exact, all of them can be systematically improved. For example, choosing N z = 16 or 32 
improves the coefficient pi to become accurate to 0.2% or 0.05%, respectively. 

The relation (|27|1 between the transverse plasma frequency and the static mass parameter 
fi 2 which determines the range of plasma instabilities is in fact exact in all cases we have 
considered, so the coefficients pk in Table HP cover both the transverse plasma frequency and 
the magnetic instability (resp. screening) mass scale for £ > (£ < 0). 



The coefficients qi j2 have been omitted out of laziness. 
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TABLE II: A comparison of parameters characterizing the dispersion laws in anisotropic plasmas 
at small anisotropy parameter £ between continuum and v-discretized hard loops. The upper part 
of the table refers to discretization by regular polyhedra M = 6, 8, 12, 20 with different orientations 
and hence different discrete (Z^) rotational symmetry around the z-axis and different numbers N z 
of distinct values v z . The numbers on the right-hand side give the factors by which the coefficients 
defined in (|A1|) . (|A2|) . (jA3|) . (|A4|) differ from the exact ones. These factors are given only up to 
the first term where there is a deviation from the exact result. The lower entries correspond to 
the "disco-ball" discretization H16|) with regular spacing in z and ip (there the results actually only 
depend on N z if Np > 2). 
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APPENDIX B: LATTICE EQUATIONS OF MOTION 

1. 1+1-dimensional case 

In the dimensionally reduced situation arising for initial conditions depending only on z 
and t, we have adjoint scalar fields <p x = A x and <p y = A y as well as auxiliary fields W v in the 
adjoint representation, which live on the sites of the spatial lattice with lattice spacing a. 
These are represented as N c x N c traceless Hermitian matrices belonging to the Lie algebra 
of SU(A / " C ). Furthermore there are the conjugate momenta ir x and n y (corresponding to d t A x 
and dtAy) and E z = —dtA z which are all represented as N c x N c traceless Hermitian matrices 
living on the temporal links. Finally, there are SU(A^ C ) group elements U which live on the 
spatial links. The latter represent the discretized version of the parallel transporter which 
in fundamental representation is given by 

U(z', z)=V exp i^-ig J* d z"A z (z")j (Bl) 

and where V represents path ordering (with z' to the left and z on the right). The discretized 
version then transports from the site z to the site z' — z + a — > s + 1, so that 



U s+ i = exp (-igaA z>s ). 



(B2) 
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The covariant derivative D z may then be discretized into a right-transporting and a left- 
transporting version, 

U], i(f> s +iU s+ i - (f) s 
-> ? , (B3) 



D 



Di 



whereas the second order derivative is then symmetric since 

jj2 _^ ^z ^z 



(B4) 

Since the E ZyS are living on the spatial links from s to s+1, the first order covariant derivative 
has to be taken as D^, while the current is replaced by an average between the beginning 
and the end of the link. Specifically, 



7T a , s (t + |) = K a , s (t - |) + e [a 1 (D L Z - Df) (j) s - g 2 [(j)a y 



a Z \ 2 /j t 



(B5) 

where a runs over x, y and a denotes y and x. The fields are then updated using symmetric 
first order spatial derivatives (D z — > D z = + D^)/2), 



4>a,s( t + € ) = 4>a,s{t) + elicit + -e) 



W v , s (t + e) = W v , s (t) + e 



—a v v y 7i y (t 



-ig[A{t), Wv(t)] + KD s z A{t) - v z D^W v>s {t) - a v v x ir x {t 



E z , s (t + f ) + U s _i(t)E z>s ^(t + |)Z7+ i(t) 



■(w 2 a v + 6v) 



•A(t) = M>^) + fA(t). (B6) 

Finally, choosing i? 2jS to transform at the spatial lattice site s rather than s + 1 we find for 
the following update law for the links U : 



U s+ i(t + e) = U s+ i(t) exp ( igeaE ZjS (t + -] 



(B7) 



Gauss's law becomes 



ijEia)> - |)] - - f ) + ^ E Wv,.(t) = (B8) 

while the approximately conserved energy (J2J) is discretized by 



£ = aE -tr 



E ZjS (t-~)+E z , s (t + ~ 



E 



1 

-tr 

4 



(vr a , s (t-^) + vr aiS (t + |)^ 



+tr (pf <M*)) 2 )] + tr <M) 2 ) 4 + e E tr \-7r XyS (t' - |) (^(f - e) + 



f=0 



&,.(t'-e)+.MO) 



(B9) 
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As initial conditions we choose 

U s+ i{0) = l Nc , ajS (O) = O, £,, s (-f) = 0, W v , s (0) = 0, 

«.(-§K*(-i)> = <T 2 5°*6 aP 6 ssl /a . (BIO) 

Under local gauge transformations A(z) one has 



4>a,s - 




= A^^Aj 1 


^a,s 




= Ag^^A^ 1 


E z>s - 




= A^^A; 1 


w v , s - 




= A.Wv^A; 1 






= A s+1 f/ s+ iA; 



(Bll) 



which can be used to check explicit gauge invariance of the results. 



2. 3+1-dimensional case 

For the 3+1 dimensional simulation, we found it convenient to follow very closely the lat- 
tice equations of motion used in Ref. [3Jj. The dynamical variables are taken to be the three 
electric fields E^x) = —d t Ai(x), lattice links Ui(x) = exp (— igaAi) and the auxiliary fields 
W v (x). Sticking to the convention of Ref. (3]]], this unfortunately creates a slight inconsis- 
tency with the definition of the links Ui in the last subsection: for the three dimensional 
simulation, these are now defined to transport from site x to site x — aii. Furthermore, we 
make the field variables dimensionless by absorbing the lattice spacing as well as the strong 
coupling constant g by 

gaA -> A, ga a E -> E, ga 3 W -> W, ga 3 j -> j. (B12) 

Making use of the shorthand notation x + aej — > x + i, the lattice equations of motion then 
read for gauge group SU(2) 

E ?(t+Lx) = Ei(t-^,x 



+ 2te tr I T a U l (t,x)J2Sj, j (t,x)) 
-« (ji(t> x ) + Ufa x)ji(t, x + i)Uj(t, x)^ 



Ui(t + e,x) = exp (i e E^t + |, x))Ui(t, x) (B13) 

W«(t + e,x) = W$(t, x) - ev ■ D s W?(t, x) + e (a v v + 6 v e,) • E a (t + ~) 

-2zefc v (u,tr (r a f/ n ^(t, x)) + V r (r a U a , yz (t, x))) , (B14) 

x) = Uj(t, x + i)U}(t, x + j)U](t, x) (B15) 



where 

is the gauge link staple, 



U a>ij (t,x) = Ui(t,x)Uj(t,x + i)U}{x + j)U}{x) (B16) 



29 



is the standard plaquette and the symmetric covariant derivative here is now 

D? (f)(t, x) = - (Ui{t, x)<f>(t, x + i)Uj{t, x) - Uj(t, x - i)<j)(t, x - i)Ui(t, x - z)) . (B17) 

Note that the sum Y^,\j\^i runs over both positive and negative directions, with U-j(t,x) = 
U}{x-j). 

Contrary to Ref. j3l| we did not discretize the first order temporal derivatives in a sym- 
metric way, thereby breaking explicit time reversal symmetry which might introduce addi- 
tional lattice discretization artifacts in the simulation. However, by using different temporal 
step sizes and different lattice spacings we have convinced ourselves that our results are 
afflicted only very mildly by these artifacts if we choose e/a sufficiently small. 

The energy density (J7|) then reads 



£ 



4X:(l-itrf/ D (t,x)) + ^tr(^(t 



e -,x) + E l (t+ t -,x)) 



2 ^' " 4^- 
e ]T toEitf - % x ) (jtf ~e,x)+ j(t', x) + Ui(t', x)j(t' -e,x + i)U}(t' , 



i t'=0 



X) 



+U i (t',x)j(t',x)U}(t',x) 



(B18) 



and Gauss's law becomes 



J2 \U}{t,x-i)E i {t--,x-i)U i {t,x-i)-E i {t 



x 



')+^E>Vv(t,x) = 0. (B19) 



Finally, we initialize our three-dimensional simulations with Ei 



2 • 



X) 



Ai{0,x) 







and 



[W«(0,x)W»(0,y)) = 5 a % y a 2 a 3 , (B20) 

where 5^ y denotes 3 Kronecker deltas, and subtract ^ v 'W a /M' 2 from each W a in order to 
obey Gauss's law. 
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